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Abstract: Hand-held OCT systems that offer physicians greater freedom to 
access imaging sites of interest could be useful for many clinical 
applications. In this study, by incorporating the theoretical speckle model 
into the decorrelation function, we have explicitly correlated the cross- 
correlation coefficient to the lateral displacement between adjacent A-scans. 
We used this model to develop and study a freehand-scanning OCT system 
capable of real-time scanning speed correction and distortion-free 
imaging — for the first time to the best our knowledge. To validate our 
model and the system, we performed a series of calibration experiments. 
Experimental results show that our method can extract lateral scanning 
distance. In addition, using the manually scanned hand-held OCT system, 
we obtained OCT images from various samples by freehand manual 
scanning, including images obtained from human in vivo. 
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1. Introduction 

Optical coherence tomography (OCT) is a high resolution optical imaging modality widely 
used in biological and medical fields [1,2]. For many clinical or intraoperative applications, a 
hand-held OCT system could be particularly useful; it would offer physicians greater freedom 
to access imaging sites of interest [3-10]. In a hand-held OCT system, it is desirable to have a 
robust and lightweight probe which can image detailed anatomical structures with a large 
field-of-view. 

In conventional OCT systems, a mechanical scanner steers the OCT probe beam to 
perform lateral scans. Sequentially acquired A-scans are assembled according to a pre-defined 
raster [1] or circumferential [10] scanning patterns to form two dimensional (2D) or three 
dimensional (3D) images. Scanners used for OCT include galvanometer-mounted mirrors, 
piezoelectric transducers (PZT) and microelectromechanical systems (MEMS). 
Galvanometers have a high linearity and accuracy; however, they are usually bulky and 
heavy, especially in the case of 3D imaging which requires two galvanometers to perform 2D 
transverse scans. PZT scanners are smaller than galvanometers and therefore are more suitable 
for hand-held probes. However, they require a high driving voltage, which is a safety concern. 
MEMS scanners are smaller but relatively expensive, and they require relatively high voltage 
[11]. 

On the other hand, OCT scans can also be performed manually, similar to manually- 
scanned ultrasound imaging systems [12,13]. A manually-scanned OCT probe without any 
mechanical scanner to steer the beam could be much simpler, cost-effective, and easy to use 
during intraoperative settings [14]. It has been shown that a simple ID, hand-held OCT probe 
integrated with standard surgical instruments can be used for 2D OCT imaging and depth 
ranging during surgery [8,15]. When surgeons manually scan the OCT probe integrated with a 
surgical tool across the target transversally, the time-varying A-scans can be acquired 
sequentially and can be used to form pseudo B-scan images. Due to the non-constant scanning 
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velocity of the surgeon's hand, such a pseudo B-scan results in a non-uniform spatial 
sampling rate in lateral dimension. Such artifact varies widely between surgeons depending on 
the stability and dexterity of their hands. 

Researchers in the ultrasound community have developed various methods in the last 
decade to correct the artifact induced by the non-constant scanning velocity in manual 
scanning, and ultrasound imaging systems have benefited from the use of manually scanned 
probes. In addition, methods including position tracking and speckle decorrelation have 
recently been adopted by the OCT community [9,14,16,17]. The speckle decorrelation 
algorithm is particularly interesting and was demonstrated a few years ago by A. Ahmad et al. 
in OCT systems for the first time [14]. Compared to a video position tracking system, the 
speckle decorrelation technique may achieve better accuracy because the dimension of OCT 
speckle is in the order of micrometers [18]; which is sufficient for high-resolution OCT with a 
micrometer-resolution. Speckle decorrelation algorithm is attractive also because it does not 
require extra hardware components and is easy to implement. 

In this study, we incorporated the theoretical speckle model into the decorrelation function 
to explicitly correlate the cross-correlation coefficient (XCC) to the lateral displacement 
between adjacent A-scans. We performed a series of experimental calibrations to validate our 
model and to show that lateral displacement between adjacent A-scans can be extracted 
quantitatively based on XCC. With the displacement extracted, we were able to correct the 
artifact induced by the non-constant scanning velocity. To test the method, we built and 
demonstrated a freehand-scanning OCT system capable of real-time scanning speed 
correction — for the first time to the best our knowledge. Our system consists of a simple 
hand-held probe, a spectral-domain OCT engine, and software for image reconstruction and 
scanning speed correction. We integrated a single-mode fiber with a needle to serve as our 
probe. The spectral domain OCT engine uses a line scan CCD camera for spectral 
interferogram acquisition. We also developed our high speed software for real-time signal 
processing based on a general-purpose graphic processing unit (GPGPU). To demonstrate the 
system, using a simple 22 gauge needle integrated with a single-mode fiber probe, we 
obtained OCT images from various samples by freehand manual scanning, including images 
obtained from human in vivo. 

One significant difference between this work and Ref [14] is that our method can directly 
calculate lateral displacement from the value of cross-correlation coefficient based on the 
speckle model we derived. Moreover, in this work, we have developed real-time scanning 
speed correction algorithm, because the scanning speed correction algorithm proposed in the 
manuscript could be easily parallelized and implemented using GPGPU. 



In this manuscript, we use a Cartesian coordinate system (x, y, z) to describe the 3D space, z 
indicates the axial direction; x is the lateral direction or the direction of the manual scan. For 
simplicity, we assume the motion of the OCT needle probe is limited to x-z plane and the 
specimen is static. 

2.1. Manually scanned OCT imaging 

As shown in Fig. 1, when a simple hand-held OCT probe without mechanical scanner is 
scanned manually in x direction, the displacement between adjacent A-scans, Ax, is a function 
of the instantaneous scanning velocity v and the A-scan acquisition rate f A , as shown in 



v varies with time for a manual-scan OCT probe; f A is usually a constant for conventional data 
acquisition devices such as a frame grabber synchronized with an internal, periodical trigger 
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2. Theory 



Eq. (1). 
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signal. As a result, Ax varies with time in the same manner as v. Therefore, the lateral intervals 
between different A-scans are different for manual scan. 



Fig. 1 . Schematic of manual-scanned OCT imaging with time-varying scanning speed. 

According to Nyquist theorem, the sampling rate, R, has to be larger than twice the highest 
spatial frequency of the specimen (F n ): R = 1/Ax = fpJv>2F n . Therefore, the scanning speed 
has to be smaller than v m shown in Eq. (2). 

v (2) 

Equation (2) also implies that a scanning velocity smaller than v m would lead to 
oversampling and information redundancy. Under the oversampling condition, there is 
correlation between adjacent A-scans. The degree of correlation can be measured by Pearson 
cross-correlation coefficient (XCC) shown as Eq. (3). 

_ ([^.y (Z) " {h.y U))] (Z + AZ) - (l x+ix , y+Ay (Z + AZ))]) 

ft,,,b).W,«,( ItAl l = ~Z ~Z ^ 

In Eq. (3), < > indicates to take the mean value of a signal. Here I xy {z) is the intensity of an 
A-scan at (x,y). I^ y (z) is calculated by taking the square of the amplitude of the A-scan. Denote 
the complex valued OCT signal as S xy (z); then I xy (z) = S x ,y(z)S* x , y (z). Similarly, h+^y+Ayiz + 
Az) is the intensity of A-scan that is displaced by (Ax,Ay,Az). <Ji xy ( Z) and £r &+Aj . V+A> . (;:+A; .) are the 
square roots of variance for I xy (z) and I x +&x,y+Ay(z + Az). 

As we assume the scanning is in x direction, Ay = Az = 0, I x +Ax, y +Ay(z + Az) becomes 
I x +Ax, y (z) and Eq. (3) becomes: 



U) " {h,y (2))] [ W (Z) " ( W (Z)}] 



%,j.(z)./» + 4,,j.(2) _ „ 



For simplicity, we use p to denote Pix, y (z),i x +Ax, y <z) m subsequent equations. 

If we assume the specimen has a homogeneous distribution of scatterers with a uniform 
scattering strength [19], e.g., the speckle is fully developed, the following relationship exist: 
<h, y (z)> = <I X+Ax , y (z) > = Iq\ <I x ,y(z) > = <I X+ Ax,y(zf> = /rms 2 - Therefore, we have: 

< ;b) =([^(z)-(^(z))] 2 ) = (/,„(z) 2 >-(/„(z)> 2 = I 2 ms -II 

<^ M = ([ (z) - (l x+M , y (z))] 2 ) = {l x+Ax , y (z) 2 ) - (7,^, (z)) 2 = 7^ s - 7 0 2 

([',,, (z) " (/„ (z)>] (z) - (J WM , (z))]) = (/„ (z)I x+Ax , y (z)> - 7 0 2 

Based on the above relationships, we can simplify Eq. (3) to: 
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P = ~ 



I 2 -I 2 
1 RMS 1 0 



Similar to ultrasound images with fully developed speckle, with the moment theorem for 
the jointly zero mean, Gaussian random variables, and assuming that the real and imaginary 
parts of S are uncorrected, we have [12,13,20], 



S X M)S\ (z) 



+ K (4) 



It is worth mentioning that to derive Eq. (4), we only utilized statistical properties of the 
random variable involved; such statistical properties are exactly the same for OCT and for 
ultrasound imaging in [12] and [20]. On the other hand, such statistical property is not related 
to the physical mechanisms of OCT image formation; and therefore the Eq. (4) is applicable 
to OCT signal. 

In the Eq. (4), hi 2 is the square of the amplitude of a complex value. Signal S X j(z) is 
determined by the physics of OCT image formation mechanism and can be expressed as the 
convolution of scattering distribution function a(x,y,z) with system's 3D point spread function 
(PSF) P(x,y,z) 

S xy (z)= J[J a[x-x\y-y',z-z')P{x',y\z')dx'dy , dz' 

x\y\z' 

Similarly, OCT signal S x+Ax<y (z) can be expressed as 

S x+Ax , y (z)= JJJ a(x + Ax-x',y-y',z-z')P(x',y\z')dx'dy'dz' 

x\y\z' 

It is worth mentioning that /indicates integration over (-oo, + oo) in the expressions of 
Sx, y (z), S x+ Ax,y(z) and in the following derivations. 

Plugging the expression of S xy (z) and S x+Axj (z) into Eq. (4) and utilizing the fact that OCT 
system's PSF is not random, we have: 



,«WO>- III III 



x,y,z x ,y ,z 



+ 1 



Assuming that the speckle is fully developed and thus scatterers in different spatial 
location are described by identical but independent random variables, we have the following 
relationship: 

(a(x-x\y-y\z-z')a(x + Ax-x",y-y",z-z")) = a 2 S(x'+Ax-x")S(y'-y")S(z'-z") 



In the above equation, a 0 is a constant representing the scattering strength. Using the 
sifting property of delta function, we have 



fff Iff 

x\y',z' x n ,y n ,z 



a 2 S(x'+ Ax- x")S(y"- y')S(z"- z') 
• P(x',y',z')P" (x",y",z")dx'dy'dz'dx"dy"dz [ 



+ 1 



J[[ a] P (x ', y ', z ') P* ( x '+ Ax, y ', z ') dx ' dy ' dz 



In OCT, the axial PSF P(z) and the lateral PSF P(x,y) are separable because axial and 
lateral PSFs are governed by different physical principles: axial PSF is determined by the 
temporal coherence of the light source while lateral PSF is determined by the imaging optics 
in the sample arm. Furthermore, in Gaussian optics model, P(x,y) is the product of PSFs in x 
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and y dimensions [18]. As a result, P(x,y,z) can be written explicitly as P(x,y,z) = P x (x) P y {y) 
P z (z) and therefore we have: 



= \al\\\p^)P y (y<)P z (z')P;{x< + te)p;(y<)p;(zyx'dy<dz\ + ll-ll 

(-co +00 +00 

\P x (x')P;(x'+Ax)dx' j P y (y')P;(y')dy' j P z (z')P z * (z')dz' 

-co J |_ —oo J |_ —oo 

2 2 2 

x-J +co +co 

= "l\p y {y')p y \y'W \ p ^')p; {z')dz' lp x ( X ')p;(x'+Ax)dx' 

—co — co —co 

Lateral PSF P x (x) can be expressed as: 



(5) 



p x ( x ) = p o ex P 



f x^ 



(6) 



In Eq. (6), w 0 is the Gaussian beam waist of probing beam [21]. It is worth mentioning that 
Gaussian beam waist in this definition is the distance from the beam axis where the intensity 
of OCT signal drops to lie. This PSF expressed in Eq. (6) is valid and consistent with 
literature [22,23]. 

Plugging Eq. (5) into Eq. (4), we have: 



\\p x (x')p;(x' + hx)ch'\ 



2 

4\ P y {y')P;{y')dy' 


2 

\p z {z')p;{z'yz' 


2 

$P x (x')P;(x'+Ax)dx' 




al\P y (y')P y '{y'yiy' 


2 +oo 

\p z {z')p;{z'yz' 


2 

j p ,( x ') p x '(x')dx' 


2 



\$p x ( X ')p;(x')dx\ 



Using the expression of P x (x) shown as Eq. (6), p can be re-written as: 



P = 



P 0 exp 



f .2 A 
X 



exp 



(x'+Ax) 



2 A 



dx 



exp 



2 x 



dx' 







( a\ 










X 






J 


P 0 2 exp 




exp 


'-4)1 






V w lj 




V w o)_ 



dx 



fexpf-2^ 



dx 



We are able to calculate the integration of Gaussian function over (-oo, + oo): 



..2 A 



[ exp -2 — j dx' = [ exp 

L [. W 0 J -oo 



2 x' 



Ax 



2 A 



dx = -^k /2w 0 



Therefore, 
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(7) 

Equation (7) shows that the value of p is merely determined by Ax for fully developed 
speckle; therefore, we can calculate the cross-correlation coefficient p between adjacent A- 
scans and use the value of p to derive the time-varying Ax as: 



Ax = w 0 lln 



1 



(8) 



It is worth mentioning that Ax with opposite signs can lead to the same value of p, 
according to Eq. (7). In other words, the scanning direction cannot be determined by 
calculating XCC. 

For digitized sample points in A-scans, p can be calculated with Eq. (9). 



Pi 



j+i 



i(v«) 2 i(w 



(9) 



In Eq. (9), i is the index of pixel in an A-scan and j is the index of A-scan. Segmentation 
of signal between if and i t is selected to calculate p. 

Although Eq. (3) and Eq. (9) implies that XCC can be either positive or negative, it is 
unlikely that XCC of two adjacent A-scans has very small or negative value when lateral 
dimension is highly over-sampled. As demonstrated previously, speckle pattern is formed by 
convolving random scattering field with OCT system's PSF. Due to the finite dimension of 
PSF, adjacent A-scans are correlated as long as lateral displacement is small compared to the 
lateral width of PSF, no matter how sample field is like. However, it is true that XCC can 
become very small or negative. To deal with this problem, in our software implementation, we 
took the absolute value of XCC for displacement estimation so that the term inside of 
logarithm operator is never negative. In addition, we applied thresholding to the value of XCC 
because small XCC indicates decorrelation between A-scans and thus is not reliable for the 
displacement assessment. It is also worth mentioning that when XCC turns negative, it does 
not represent the change of scanning direction. 

2.2. Scanning speed correction based on speckle decorrelation 

Equation (8) indicates that lateral interval between A-scans can be extracted from the XCC 
and therefore we could use XCC to correct artifact induced by non-constant scanning speed. 
The flow chart of the scanning speed correction for one frame of data is shown in Fig. 2. In 
our OCT system, the interferometric spectra are detected by a line scan CCD camera and the 
data is transferred to the host computer through a frame grabber as frames. One frame consists 
of N spectra acquired at lateral locations: Although Ax„ the interval between X\ 

and Xi + 1, is not a constant for different A-scans due to non-constant scanning speed, we could 
extract Ax, using p h the XCC between 7,(z) and I i+ \{z) (A-scans obtained at spatial coordinate 
X\ and X; + t ). As a result, we were able to estimate Ax toto/ , the displacement (in x direction) 
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between the first and the last A-scan in a frame by summing up Ax,: Ax tolal = ZAx,. With a pre- 
set sampling interval Ax s , we could calculate the number of A-scans required for this 
particular frame of data by dividing Ax tota i with Ax,: M = Ax (oto/ /Ax,. Afterwards, we performed 
interpolation to obtain A-scan data at spatial points 0, Ax s , 2Ax s ,...,MAx s to obtain A-scans that 
were evenly distributed in x dimension. 



i=l 




i<N 



calculate p t 



A. 



calculate Av, 



i=m 



1 



X 



Setup interpolation 
locations:0\ Ax s , 

I 



Interpolation to 
redistribute 
A-scans 



Fig. 2. Flow chart for the scanning speed correction using cross-correlation coefficient. 



3. OCT system and software implementation 

The basic properties of the OCT system used for our calibration experiments have been 
reported in our previous work [24]. Briefly, it was a spectral domain OCT system based on 
GPGPU processing with a 70kHz A-scan rate. In that set-up, an achromatic doublet lens at the 
sample arm was used to focus the incident light beam and collect back-scattered photons. We 
used the focal length of imaging objective, focal length of the collimator, and mode field 
diameter of the fiber to calculate wq. From our calculation, wq equaled 6(a.m and this indicated 
a 12um lie lateral resolution of our OCT. This was further verified by using our OCT system 
to image 1951 USAF resolution target. The obtained en face OCT image clearly showed that 
our OCT system can resolve the 5th element in group 6 which corresponds to a lOitm FWHM 
lateral resolution and therefore 12(j.m lie lateral resolution. Therefore, we assumed vv 0 to be 
6(a.m for this system when using Eq. (7) to correlate p and Ax. With the 12|im lateral 
resolution, F n = (l/12)((o.m _1 ); therefore the largest scanning speed that satisfies the 
requirement of Nyquist sampling is about 420mm/s, as implied by Eq. (2). For our calibration 
experiments, a high-speed galvanometer was used to perform lateral scanning. 

Our freehand-scanning OCT was a modified version of this system. As shown in Fig. 3(a), 
we used a superluminescent diode centered at 840nm with a full width half maximum 
bandwidth of 55nm (Superlum, S840-B-I-20) as a broadband source. The interferometric 
signals are dispersed spectrally by a spectrometer that consists of a collimator, diffraction 
grating (Wasatch Photonics, HD 1800 1/mm @ 840nm), and two identical achromatic doublet 
lenses (AC508-250-B, f = 250mm) for focusing. The spectra are detected by a line scan CCD 
camera (e2v, AVIIVA EM4, maximum line rate 70kHz). Spectral data is transferred through a 
high-speed frame grabber (Matrox Solios eV-CLF) into our host computer (Dell Precision 
R5500 Rack Workstation, 6-Core Intel® Xeon® Processor X5690 3.46GHz, 12GB RAM, 
64bit Windows 7 operating system). For high-speed signal processing, we implemented 
massive computation on a GPGPU (Nvidia GeForce GTX 480) with 480 cores, with each core 
operating at 1.4 GHz. To build a simple, lightweight, and small probe which can have 
arbitrary length, we adopted a common path (CP) configuration for our interferometer [25,26]. 
In the CP interferometer, the reference and sample light shares the same probe arm which is 
simply a single mode fiber with its tip cleaved in right angle, as in Fig. 3(b). The reference 
light comes from the Fresnel reflection at the fiber tip. The sample and reference light is 
routed by a 50/50 fiber optic coupler to the spectrometer. To protect the fragile fiber tip, we 
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integrated the fiber probe with a 22 gauge needle attached to a syringe so that it can be easily 
hand-held, shown as Fig. 3(c). The system performance was characterized in our previous 
work [26]. Despite the fact that Gaussian beam diverges, the lateral resolution of our bare 
fiber probe is higher than 25(j,m from results obtained from both experimental measurement 
and ZEMAX simulation [27]. 




Fig. 3. (a) OCT system diagram; (b) principle of common path interferometer based on single 
mode fiber; (c) single mode fiber probe integrated with needle. 

The signal processing procedure of our system is briefly summarized in Fig. 4. A data 
frame containing N spectra was acquired and transferred to GPU memory. Afterwards, we re- 
sampled the spectral data from wavelength (1) space to wavenumber (A:) space using cubic 
spline interpolation and then performed fast Fourier transformation (FFT) to obtain A-scans. 
With the obtained A-scans, we calculated the XCC between adjacent A-scans using the OCT 
signal intensity and re-distributed the A-scans using algorithm shown in Fig. 2 to achieve 
uniform spatial sampling. Before calculating p, we processed each A-scan with a moving 
average filter that averages three adjacent pixels in axial direction and subtracted the output of 
the filter from each A-scan to reduce low spatial frequency components in the A-scan and 
therefore increase the sensitivity of cross correlation calculation for lateral motion estimate 
[14]. Due to the high A-scan rate of our OCT system, lateral dimension was highly over- 
sampled during the manual scan; therefore simple nearest neighbor interpolation could 
achieve satisfactory result in re-distributing A-scans and achieve uniform lateral sampling. To 
match the dynamic range of the OCT data and display device, we applied a truncated log 
transform to the OCT signal before transferring the signal back to the host computer for 
display. Procedures shown as red blocks in Fig. 4 were implemented with GPU. Moreover, we 
implemented multi-thread programming so that data acquisition was in parallel with 
processing; therefore, we were able to acquire spectral data detected by the CCD continuously 
as long as we kept the data acquisition rate slightly lower than processing rate. Our software is 
able to process over 62,000 A-scans every second. Although slightly less than the maximum 
data acquisition rate of our camera (70k line rate), this speed allows a maximum lateral 
scanning speed to be approximately 620mm/s assuming wo is about 20um for the single mode 
fiber probe and this lateral scanning speed is significantly larger than moderate manual 
scanning speed (several millimeters per second). With the software optimization, we will be 
able to further increase the processing speed in future studies. 

To obtain an accurate value of w 0 in our software for our single mode fiber probe, we have 
first used an estimation of w 0 based on the experimentally measured lateral resolution of our 
CP OCT system from our previous work [26,27] and indicate this value as w. Afterwards, we 
manually scanned a highly scattering phantom for 1cm (L = 1cm) and acquired a certain 
number of A-scans that were uniformly distributed. We performed such scanning for 10 times 
and calculated the average A-scan number in all the measurements which is indicated by M. 
According to Eq. (8), we were able to obtain a better estimation of vv 0 that equals 
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(wL)/(MAx s ). We varied the value of Ax, and obtained a consistent value of w 0 . All images in 
section 5.2 were acquired based on w 0 = 18.5um which was constant for different values of 

Ax s , 



data ^^^^^^^^^^^^^^^^^^M Transfer data to 

transfer data to 
memory 



I t 




Fig. 4. Block diagram of signal processing for the OCT system with real-time scanning speed 
variance correction. 

4. Calibration of the relationship between cross-correlation coefficient and lateral 
displacement 

In our calibration experiments, we scanned a phantom consisting of 9 layers of cellophane 
tape to verify the relationship between displacement Ax and p (XCC), shown as Eq. (7) and 
(8). We used a galvanometer to perform lateral scans with known scanning speeds. We 
applied a periodical sawtooth voltage V from a function generator to the galvanometer and 
synchronized V with the acquisition of a frame of data which contained N A-scans (N = 1000). 
For a 100% duty cycle sawtooth driving voltage, Ax, lateral interval between adjacent A-scans 
stays constant because the driving voltage increases linearly during signal acquisition. 
Therefore, we could calculate the displacement between adjacent A-scans directly from the 
amplitude of the sawtooth function: Ax = yVIN. Here y is a coefficient that relates the driving 
voltage (V) applied to galvanometer and the probing beam displacement (D) at the focal plane 
of the imaging lens: y = D/V. y was measured to be 1.925mm/V in the OCT setup for our 
calibration experiments. As a result, by applying different V, we could achieve different 
scanning speeds and thus different Ax. We acquired B-scans at various scanning speeds. One 
example of the image obtained is shown 5(a) which contain 1000 A-scans, with a sampling 
interval Ax equal to 0.96um. 

For B-scan acquired at a certain spatial sampling interval determined by the driving 
voltage, we calculated p„ XCC between the i and (z + l) th A-scan, using pixels within the 
range as indicated by the double-head red arrow in Fig. 5(a). Afterwards, with all the 
obtained, we took the ensemble average of XCC: 




Using the above processing, we obtained a value of p from each B-scan corresponding to a 
certain At. The result is shown as red circles in Fig. 5(b). As a comparison, we also plotted the 
theoretical relationship between p and Ax shown in Eq. (7) as a black, dashed curve in 
Fig. 5(b). By calculating ensemble average of XCC using A-scans with different offsets from 
the same B-scan, we obtained decorrelation curve as shown in Fig. 5(c). Similarly, theoretical 
relationship between p and Ax is shown as black, solid curve in Fig. 5(c). The consistency 
between experimental results and the analytical model described by Eq. (7) and (8) implies 
that we can use p to quantitatively extract the lateral sampling interval and thus correct non- 
constant scanning speed. 

XCC, as the correlation of two random variables, is inherently a random variable. It is 
very important to evaluate the statistics of XCC to assess the accuracy in displacement 
estimation. To evaluate the statistics of XCC, we calculated the standard deviation and the 
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mean of p, from B-scans acquired experimentally at different spatial sampling intervals (<&,). 
The results are shown in Fig. 5(d). With each obtained p„ we used Eq. (8) to calculate a 
corresponding displacement value Ax, and then assess the variation of Ax„ cr xl 2 . Assume that 
the displacement between each pair of the A-scans follows the same statistics and the probe 
travels for a given distance Ax tot ai- In this case, M, the number of A-scans acquired is 
Axtotai/Sxi. Based on these assumptions, a xtola i 2 , the variance of the estimated displacement 
approximately equals Ma x f. In Fig. 5(e), we show the ratio between a^ xtotal and Ax IMaI , for 
different values of Ax lota i and Sx h 

As shown in Fig. 5(b) and especially Fig. 5(c), when interval between A-scans is small 
(<5um), the measured XCC values are highly consistent with results from theoretical 
calculation. However, with larger interval (>5u.m) between A-scans, the measured XCC 
values become slightly different from the values calculated using Eq. (7). Moreover, Fig. 5(d) 
and 5(e) show that a smaller sampling interval (smaller dxi) and larger lateral distance 
travelled would result in a higher accuracy of displacement estimation. In other words, a large 
number of sampling points for a given displacement can provide a better displacement 
estimation, due to the inherent statistics of XCC. Other than the random nature of XCC, there 
are several reasons for displacement calculated using XCC and actual displacement to be 
different. First, OCT signal suffers from optical and electrical noise, such as shot noise, excess 
noise, thermal noise, etc; in addition, OCT signal decorrelates overtime due to random 
environmental disturbance. Second, the waist of probing beam varies at different lateral 
displacement due to lens abbreviation. Third, parts of A-scans with low signal intensity 
produced high correlation due to signal absence, so that the measured XCC was higher than 
the theoretical values for large displacements in Fig. 5(b) and 5(c). However, those factors 
might be negligible as compared to the inherent random noise of XCC, because OCT is a high 
speed and high sensitivity imaging modality. 

Results in Fig. 5(d) and 5(e) implies that errors in displacement estimation are minimized 
with small sampling interval which can vary over time and oversampling in x dimension is 
necessary for accurate scanning speed correction. Comprehensive evaluation of the error in 
distance measurement using XCC will be our future work and is beyond of the scope of this 
manuscript. 




Sx t ( fim ) 5x t ( fim ) 



Fig. 5. (a) Image obtained in calibration experiment (Ax = 0.96um); (b) relationship between p 
and Ax obtained by calculating XCC between adjacent A-scans from different B-scans (red 
circles: experimental; black, dashed line: theoretical); (c) relationship between p and Ax 
obtained by calculating XCC using A-scans with different offsets from the same B-scan (red, 
solid line: experimental; black, dashed line: theoretical); (d) the ratio between standard 
deviation and mean of p t , at different sampling intervals; (e) ratio between CTawhoI and Axtoai. 
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5. Results 



5.7 Quantitative lateral sampling interval extraction 

We used the same OCT setup for this experiment as in the calibration experiments. To 
demonstrate that we can use XCC for the quantitative displacement extraction and thus 
quantitatively correct artifacts from non-uniform scanning speed, we applied a sinusoidal 
driving voltage (f = 28Hz, Vpp = 1.5V) to the galvanometer and scanned the multilayered tape 
phantom. Knowing voltage applied to galvanometer and the value of y from our previous 
measurement, we were able to calculate the displacement of the probe beam with regards to 
its neutral position when the voltage applied to the galvanometer was 0. The calculated 
displacement of the probe beam at different time is shown in the upper inset of Fig. 6(a). 
Instantaneous Ax was calculated by taking the absolute value of the difference between the 
displacements of adjacent sampling points, as shown in the lower inset of Fig. 6(a). In this 
experiment, 5000 A-scans were acquired sequentially. By stacking the A-scans, we obtained 
the pseudo B-scan shown in Fig. 6(b). When the driving voltage/displacement reaches their 
extreme points, the interval between adjacent A-scans became smallest, which can be clearly 
seen in Fig. 6(a). With a small sampling interval, data is redundant, as in areas indicated by 
the arrows in Fig. 6(b). We calculated the XCC between the adjacent A-scans in Fig. 6(b). 
XCC as a function of time is shown in the upper inset of Fig. 6(c) which was processed with a 
low pass filter for noise reduction and normalized to the maximum value. We further 
calculated Ax using Eq. (8) with XCC obtained in the upper inset of Fig. 6(c) and show the 
result as the red curve in the lower inset of Fig. 6(c). In the lower inset of Fig. 6(c), we also 
plot Ax calculated from the known driving voltage applied to the galvanometer, as the black 
curve. The consistency between the red and black curve verified our assumption that Ax could 
be extracted from XCC quantitatively. There are several reasons for the red and black curve in 
the lower inset of Fig. 6(c) to be slightly different, as discussed in the previous section. 
However, the inherence statistics of XCC plays the most significant role in resulting errors. 
As shown previously in Fig. 5(d), if the sample is laterally homogenous and the displacement 
between the adjacent A-scans is small, for example less than lum, the errors in displacement 
due to the inherent randomness of XCC are small. Therefore, the errors might be mostly due 
to other random noise in OCT signal. With a larger sampling interval, errors come from the 
inherent statistics of XCC rather than other noise in OCT signal. As shown in the lower inset 
of Fig. 6(c), difference between estimated and actual interval is smaller when interval between 
A-scans is smaller. 

To validate the scanning speed correction algorithm, we took A-scans between 8ms and 
13ms when the scanning velocity did not change its direction; thus, we didn't have to consider 
the ambiguity of scanning direction. Simply stacking A-scans acquired, we obtained Fig. 6(d) 
which exhibits an oversampling artifact on the left side of the image, as indicated by the 
arrow. To remove such artifact, we set sample interval Ax, to be 5|im and performed nearest 
neighbor interpolation as described in section 2.2. The resultant image is shown in Fig. 6(e) in 
which the oversampling artifact is removed. 
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Fig. 6. (a) Displacement of probing beam versus time (upper); Ax as a function of time (lower); 
(b) pseudo B-scan obtained from sinusoidal scanning pattern; (c) upper inset: XCC calculated 
from adjacent A-scans; lower inset: Ax calculated from XCC (red) and ground truth Ax 
calculated from the driving voltage (black); (d) pseudo B-scan with artifact induced by non- 
constant scanning speed; (e) B-scan after non-constant scanning speed correction. 

To demonstrate the effectiveness of our method more clearly, we manually scanned a 
phantom consisting of multiple layers of tape and saved all the A-scans obtained. Stacking all 
the A-scans, we obtained Fig. 7(a) which suffers from motion artifacts, as indicated by red 
arrows. After correcting A-scans using the method proposed in this paper, we obtained 
Fig. 7(b) which is free of distortion due to non-constant scanning speed. It is worth 
mentioning that scale bars in Fig. 7 are only applicable to axial dimension. 
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Fig. 7. OCT images obtained from manual scan: (a) before scanning speed correction; (b) after 
scanning speed correction. Red arrows in Fig. 7(a) indicates areas with motion artifacts. 

5.2. Images obtained from manually scanned OCT probe with real-time correction 

Using our scanning speed corrected, hand-held OCT system, we manually scanned our single 
mode fiber probe across the surface of an infrared (IR) viewing card by moving the probe 
with a freehand. In our real-time scanning speed correction software, we set the spatial 
sampling interval Ax, to be lum, 2um, and 4(j.m and show the corresponding images obtained 
in Fig. 8(a), 8(b), and 8(c). The plastic covering film and the underneath fluorescence 
materials of the IR card can be clearly seen from Fig. 8. With different spatial sampling 
interval Ax s , the same physical length is represented by different numbers of A-scans. As the 
lateral axis of Fig. 8 is A-scan index, the scale of porous structure of the fluorescence 
materials decreases from Fig. 8(a) to 8(c) due to the increasing sampling interval. Results in 
Fig. 8 verify that we were able to achieve uniform spatial sampling interval during manual 
scan and the sampling interval is explicitly determined through Ax, a which is parameter in 
our software. 




Ascan index Ascan index Ascan index 



Fig. 8. Images of IR viewing card with sampling interval Ax, equal to lum (a), 2um (b), and 
4um (c). 

To further evaluate the accuracy of our scanning speed correction method, we imaged a 
quality resolution chart with 1 line per mm from Edmund Optics through manual scan, as 
shown in Fig. 9(a). The red arrow in Fig. 9(a) indicates the scanning direction. Figure 9(b) is 
the image obtained from the software with real-time correction capability. Periodical structure 
is clearly shown in Fig. 9(b). To quantitatively evaluate the accuracy of our re-sampling 
algorithm, we averaged signal amplitude of each A-scan in Fig. 9(b), performed mean 
subtraction and obtained the blue curve (M,) shown in Fig. 9(c). Afterwards, we extracted 
zero-crossing points of M- l to detect the edge of each line, as indicated by red circles in 
Fig. 9(c). Then we calculated widths of the lines, their mean and standard deviation (STD). 
The ratio between STD and mean of the width was 0.025, which indicates that our method 
effectively removed artifacts induced by non-constant manual scanning speed. In comparison, 
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in Fig. 9(d), we show the image obtained from manual scan, but without correction using 
cross-correlation. Artifacts due to non-uniform scanning speed are clearly visible in Fig. 9(d). 



BO Edmund 



lity Resolution Chart 





Fig. 9. (a) photo of quality resolution chart; (b)OCT image obtained from manual scan with 
scanning speed correction; (c) Blue curve: Mean OCT signal of difference A-scans in Fig. 9(a); 
red circles: zero-crossing points of the blue curve; (d) OCT image obtained from manual scan 
without scanning speed correction 

We have also manually scanned the skin of a healthy volunteer using our hand-held OCT 
probe with 2u.m digital sampling interval. To perform manual scan, one of the author held the 
probe almost perpendicular to the sample surface and moved the probe laterally. Images 
obtained from the index finger tip and the palm are shown in Fig. 10(a) and 10(b), 
respectively. White scale bars in Fig. 10 represent lOOum and arrows indicate sweat duct. 
Epidermis and dermis layer can be visualized in Fig. 10. As light can penetrate deeper in the 
palm skin, the subcutaneous layer can also be seen in Fig. 10(b). 
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Fig. 10. Manually scanned OCT image of human skin from finger tip (a) and palm (b). 

To further demonstrate the effectiveness of our method on heterogeneous samples, we 
performed manual scan on an onion sample and show the obtained image as Fig. 1 1 in which 
hexagon shaped onion cells can be observed. 
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Fig. 1 1 . Image of onion cells obtained from manual scanning. 



6. Discussion 

Equation (7) forms the mathematical foundation of this work. In deriving Eq. (7), it assumes 
that the speckle is fully developed. Therefore, to validate Eq. (7) experimentally, we used 
several different models to test our method. We used a multi-layered phantom without much 
heterogeneity in lateral dimension for our calibration experiments. However, most OCT 
images using real specimens have partially developed speckle instead of fully developed 
speckle. If the sample is heterogeneous, the correlation coefficient between adjacent A-scans 
not only depends on the lateral interval, but also depends on sample structure. Moreover, 
when the probe scans across a boundary within the sample, due to the abrupt change in OCT 
signal, a new A-scan will be attached to the data set regardless of the lateral displacement 
between the current and previous A-scans. As a result, heterogeneous sample can cause 
inaccuracy in lateral motion correction of our method. However, for highly scattering samples 
such as skin, it is usually reasonable to assume that areas corresponding to sample boundary 
take only a few pixels and therefore do not contribute significantly in the calculation of XCC. 
As a result, Eq. (7) is valid for highly scattering specimens when a significant portion of the 
specimen contains homogeneous scatterers, although speckle does not develop fully in most 
biological specimens. This was verified in the experiment using a quality resolution chart with 
abrupt changes in lateral features as a sample. We further tested and verified the method using 
in situ tissue imaging. 

As indicated by Eq. (7), to quantitatively estimate Ax from XCC, it requires that we know 
vv 0 , the Gaussian beam waist of probing beam which can be calculated or experimentally 
measured. As a result, the calibrating decorrelation curve shown in Fig. 5(b) is only valid for 
an OCT system with a certain vv 0 . If w 0 used in the calculation for lateral interval between 
adjacent A-scans is different from the actual beam size, the image reconstructed from our 
algorithm will be different from the "true" image by a scaling factor in the lateral dimension. 
However, under such circumstance, uniform sampling can still be achieved to obtain an image 
that is easy for human to comprehend. In fact, the size of the imaging beam from the probe 
changes as the beam propagates, and the lateral PSF of OCT system depends on the image 
depth. Therefore, the speckle decorrelation curve has depth dependency as well. Considering 
the overall effect, the lateral resolution defined by the Gaussian beam waist is always slightly 
different from the decorrelation length of OCT signal. Moreover, to reduce the effect of the 
diverging beam, we took only part of an A-scan to calculate XCC when implementing our 
software. As a result, the statistics of different pixels within the segment of an A-scan does 
not vary significantly. 

In this work, our assumption is that the manual scan is one dimensional in x axis and that 
there is no axial motion from the scanning probe, which is not exactly true in a realistic 
scenario. For example, human hands suffer from physiological tremor and this makes the 
probe to move randomly in both lateral and axial directions. However, experimental results 
have shown that the tremor during retinal microsurgery has low frequency motion (<lHz) 
with amplitude in the order of lOOum [28]. As a result, with high data acquisition rate, 
adjacent A-scans usually do not have offset in axial direction for more than a few pixels. 
Moreover, a cross-correlation maximization-based shift correction algorithm was recently 
proposed to suppress artifact due to axial motion [29], which might be helpful to improve the 
performance of our image acquisition algorithm in the future. In our future study, we are 
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going to conduct a more comprehensive theoretical and experimental study on motion 
tracking using OCT speckle texture analysis, by considering different types of motion 
including axial translation, lateral translation and rotation. 

Another drawback of our method is that it is not able to determine the direction of 
scanning; therefore a correct reconstruction of sample structure requires that the scan is in a 
single direction. This might results in ambiguity if the surgeon's hand moves back and forth 
around a region of interest. The problem can be solved by combining OCT scan with video 
microscope technology. The spatial location of OCT tool can be tracked in videos and 
therefore the scanning direction can be obtained with slightly lower time accuracy. 
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